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Buneman instability is often driven in magnetic reconnection. Understanding how 
velocity shear in the beams driving the Buneman instability affects the growth and 
saturation of waves is relevant to turbulence, heating, and diffusion in magnetic 
reconnection. Using a Mathieu-equation analysis for weak cosine velocity shear to- 
gether with Vlasov simulations, the effects of shear on the kinetic Buneman insta- 
bility are studied in a plasma consisting of strongly magnetized electrons and cold 
unmagnetized ions. In the linearly unstable phase, shear enhances the coupling be- 
tween oblique waves and the sheared electron beam, resulting in a wider range of 
unstable eigenmodes with common lower growth rates. The wave couplings gener- 
ate new features of the electric fields in space, which can persist into the nonlinear 
phase when electron holes form. Lower hybrid instabilities simultaneously occur at 
K/k± ~ ■>/ m e /rrii with a much lower growth rate, and are not affected by the velocity 
shear. 
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I. INTRODUCTION 



The Buneman instability has been extensively studied in theory and simulation since it 
was discovered in 1958^. It is well-known in one- dimensional (ID) theory that once the 
relative drift of ions and electrons exceeds the threshold of approximately twice the electron 
thermal velocity, the interactions between waves and electrons will lead to the growth of 
the Buneman instability 2 -^. Lower hybrid instabilities lie on an oblique branch of current- 
driven instabilities^ with k„/k ~ v/m e /mj. Their relation to Buneman instabilities has not 
been clarified yet in two-dimensional (2D) spectral space. Recently, interest in Buneman 
and lower hybrid instabilities has been renewed because of their importance in magnetic 
reconnection^— . Magnetic reconnection is one of the most relevant mechanisms associated 
with explosive events in nature and in laboratory experiments, such as solar flares, substorms 
in the magnetosphere, and sawtooth crashes in fusion experiments. 

Magnetic reconnection can convert magnetic field energy into thermal and kinetic energy. 
Oppositely directed magnetic fields merge and lead to the release of a significant fraction of 
the stored magnetic energy, which produces two regions of fast outflow (see Fig. [T|). One of 
the most important problem in understanding magnetic reconnection is determining what 
mechanisms can facilitate reconnection fast enough to explain the explosive events observed. 
From magnetohydrodynamic (MHD) models, the convection and diffusion of the magnetic 
field can be described in terms of the magnetic Reynolds number R m = 47rL Vo/c 2 ?7 2 , where 
Lq and Vo are respectively the typical plasma velocity and magnetic field spatial length. 
If R m <C 1, the collisional resistivity can effectively dissipate the magnetic field's energy 
and facilitate fast reconnection. However, the following question remains: What drives fast 
reconnection if the collision-induced resistivity is not sufficiently large (i.e. R m > 1)? This is 
a condition common in both space and laboratory plasmas. One of the most promising and 
physically interesting mechanism proposed to answer this question is turbulence-induced 
dissipation facilitating fast reconnection^— . 

For turbulence- induced dissipation to work, the turbulence has to be generated sponta- 
neously during the magnetic reconnection. In the diffusion region of magnetic reconnection, 
current sheets form as a result of changes in the magnetic field configuration. The intense 
thin current sheets which develop in strong guide field reconnection can drive streaming 
instabilities&ii^kL^ and electron velocity shear instabilities 2 ^— etc. The latest 3D PIC 
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simulations show that anomalous momentum transport generated by an electromagnetic 
electron velocity shear instability can influence the magnetic reconnection process^, but if 
turbulence resistivity can affect reconnection rates is still unclear. A deeper understanding 
of current-driven electrostatic instabilities is required. Buneman and lower hybrid instabil- 
ities are two of the most common electrostatic instabilities driven within current sheets in 
magnetic reconnection. They can produce electron holes and strong turbulent heating near 
the x-line and near the separatrix— >2^rM. (Fig{T]). 

Current sheets become thinner and thinner during reconnection, and velocity shear is 
generated in the current sheets regardless of the initial velocity distribution. Extensive 
studies on new instabilities driven by electron velocity shear within current sheets have been 
performed^— i2£isE However, only a few studies have been carried out regarding how velocity 
shear affects the classical Buneman and lower hybrid instabilities 33 although their role is 
important to the understanding of turbulence-induced dissipation processes in reconnection. 

In this paper we investigate the role of weak velocity shear on Buneman and lower hybrid 
instabilities in 2D. Specifically, we study a plasma model consisting of strongly magne- 
tized electrons and cold unmagnetized ions. A Mathieu equation analysis, first proposed by 
Goldman 3 ^ 1 ^, is used and the results are compared with those of Vlasov simulations 3 ^ 1 ^. The 
results we obtain from the Mathieu equation analysis are found to be consistent with the 
results from the Vlasov simulations with weak initial velocity shear, validating the approxi- 
mations we adopted in the analytical method. We compare the Buneman instability driven 
by electrons with a uniform velocity drift (i.e., the classical Buneman instability) with the 
instability driven by drifting electrons with a weak cosine velocity shear. We find that the 
Buneman instability is no longer a purely magnetic field-aligned instability as commonly 
assumed in uniform beam^" 3 ^. The velocity shear enhances the interplay between oblique 
waves and electrons and produces a wide range of eigenmodes with a common growth rate 
lower than in the uniform-drift case. As shear increases, the fastest growing mode changes 
from a parallel plane wave to an eigenmode with significant oblique wavenumber content, 
as discussed in Appendix [A] The shear does not significantly change the growth rate of 
the co-existing lower hybrid instability, which is much weaker than that of the Buneman 
instability under our assumptions. We obtain the eigenfunctions in 2D spectral space from 
the Mathieu equation analysis. These eigenfunctions show that the shear greatly modifies 
the 2D spatial structures of the instability-induced electric fields. 
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FIG. 1. Illustration of magnetic reconnection. The red color indicates the electron diffusion region 
around the x-line with a scale of the electron skin depth c/uo pe . The semi-opaque blue color indicates 
the ion diffusion region with a scale of the ion skin depth c/u p i. The inflow velocity is Vi and the 
outflow velocity is v a . The stars "*" mark some of the possible regions where electron holes can 
develop. 

II. BASIC EQUATIONS AND SOLUTIONS 

A. Basic equations 

We assume 1) the electrons are strongly magnetized such that r e = v te /fl e "C A e = v te /cu pe , 
where r e is the electron cyclotron radius, A e is the Debye length, u pe is electron plasma 
frequency, v te = \JkT /m e is the electron thermal velocity, and k is Boltzmann's constant. 
2) the kinetic electrons are constrained to move along magnetic field lines while the ions are 
treated as unmagnetized and cold (i.e., the ratio of ion and electron temperature satisfies 
Tj <C T e ). The validity of this assumption requires that the wave phase speeds are larger 
than ion thermal velocity. 

We study the instabilities in the 2D z-y plane. The electrons move along the z direction 
with drift velocity Vb(y). The magnetic field is treated as infinite in the z direction for 
electrons and zero for ions, with no initial electric field (Eq = 0). The initial ion drift 
velocity is also zero. The electron drift velocity is modulated by a weak cosine velocity 

4 



shear: 



v b (y) = v [l + ecos2n (y/L y - 0.5)}, y e [0,L„], (1) 

where e is a small quantity and L y is the (spatial) periodicity of the shear profile. 

Because the velocity shear is a function of y, the perturbed electric fields generated 
by instabilities are also functions of y. With the assumptions of cold ions and infinitely 
magnetized electrons, the perturbed and unperturbed functions do not depend on v y . We 
perturb the Vlasov equation with f q = f 0q + Sf q = f 0q (v qz ,y) + 5f q (v qz ,y)e l( - k * z ~ UJt \ where q 
represents electrons (e) or ions (i). The electrostatic potential is <5$ = <fi(y)e' l( - kzZ ~ u;t \ The 
cold (|w/fc| > v t f, here v t i is the ion thermal velocity) initial ion distribution function can be 
written as fa = noS(vi y )S(vi Z ). To properly build our model, the initial electron distribution 
function is approximated by a ID drifting kappa-function 39 with k — 1: 



JOe 



It o 



J te 



(2) 



The kappa-functions used to represent the electron distribution are quasi-Maxwellian at 
sub-thermal velocities, but with power-law rather than exponentially decreasing tails at 
supra-thermal velocities. Such power-law tails are a common feature of measured distribu- 
tions in collisionless space plasmas. In our model the kappa distribution also simplifies the 
functional form of the susceptibilities, and thus enables us to obtain a specific mathematical 
differential equation with general solutions. We have compared our theoretical predictions 
to the results of 2D Vlasov simulations using both Maxwellian and kappa electron distri- 
butions, and find that the behavior in the two simulations is qualitatively equivalent. Only 
results from the simulations with Maxwellian electron distributions are presented in this 
paper. 

After inserting the perturbations and initial conditions into the first order linear Vlasov 
equation and Poisson equation, 

ddfg ^ d6fg _ ^dS^df^ = Q ,^ 
dt qz dz m dz dv qz 

- A5<5>(y,k z ,u) = 4vr(p e + Pi ), (4) 
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we obtain 

& 



W pe f f 0e (v ez )dv 

ez / n \ 



Upj f f i(viz)d v, 

where C indicates integration along the Landau contour in the complex plane, and p e and 
Pi are the electron and ion charge density respectively. 
Upon substituting f 0e and fa into ([6]) and ([7]) we find 

I ! 2 

_ _^pi_ / Q \ 

Xizz n > \o) 

Xezz fc^tf+o 3 ' lJj 

where f = [w - 

We normalize the quantities by defining w — )■ a) = u/u pe ; k z K z = k z vo/u pe ; v tez — > 
u e = vtez/v ; S = {L y u pe / '-kv ) 2 , so that 

^_ SA1(1 + 0, (10, 
Xi 22 = -wj/w 2 , (11) 

where 6* = it(y/L y — 0.5), G [— it/2, it/ 2], For later use, we also define the normalization 

ky Ky = kyV Q /UJ pe . 

Finally, we expand \ ezz in e to the first-order and obtain 

_| + ( a _2g cos 2^)0 = 0, (13) 
where is a function of 9, the parameters a and q are complex, and defined as 

V = ^^, a = -SKl{l + -^) (14) 

-fr 2 «e 1 + Xizz 



and 



so that 



XezzO Xezz\£=0 T ^ / .no- K*-®) 

Kjuj, (r) + iy 



Equation f JTB"]) is the well-known Mathieu equation. 



6 



B. Solutions of the Mathieu Equation 



The Mathieu equation f lT3|) is written in standard forrn^. The quantity a is the char- 
acteristic value (or eigenvalue) and the parameter q is defined in (fT5|) . For specific pairs 
(a, q), the Mathieu equation has a unique analytical solution (eigenf unction) <p, which can 
be written as <fi r {9) = e ird <j) r (9) , where r is an integer or a rational number. The value r is 
called the characteristic exponent, and r is a complex function of 9 that can be either even 
or odd. For periodic boundary conditions, (f) r has period tt, and r is required to be integer. 
In the case of the electric field, r is more restrictively required to be an even integer. If 
q = 0, Eq. f fT3|) reduces to the standard oscillator equation, where r = ^fa and the solutions 
reduce to (fi r oc cos(y / a6 l ) (even function) and sin( A /a6 ) ) (odd function). 

The fact that <p r = e ire cj) T {9) suggests that the parameter r is related to to the perpendic- 
ular wavenumber k y . For the case q = 0, a map r = r(k y ) can be established. In equation 
({TBI , r is related to 9 = 7i(y/L y — 0.5) so that 9 scales with the width of the box L y ; thus 
k y = rix/Ly. However, for the case q ^ 0, the correspondence between r and k y is nontrivial 
and needs an alternative treatment (see Appendix [A]) . 

We solve equation (|T3l) with the given periodic boundary conditions in r-K z space for 
e = and 0.2. The case of e = implies q = 0, which corresponds to uniform velocity drift. 
For specificity, we choose the following parameter values: v = 5v te , L y = 256 A e , and the 
ratio of ion to electron mass is 1836, so that S = 266. These parameters are the same as those 
used in the Vlasov simulations described in the next section. To allow comparison between 
the analytical solutions and our Vlasov simulations, we use the electric field eigenmodes E e z 19 
and E y 9 instead of $: 

53> r = ^y^ 2 -^, (17) 
Ef 9 = -Re {^^j = -Re {ik^y^-^) , (18) 

^ = -ifc(^)=-ifc(^^), d9) 
Ey ig = ^-Ef 9 . (20) 

Ly 

For periodic boundary conditions we require not only that the eigenfunction be peri- 
odic, but also require that Ey 19 be periodic so that the electric fields are continuous at the 
boundary. As we have mentioned only the eigenfunctions with even integer r can satisfy 




FIG. 2. Eigenmodes of E z and E y for z = and t = calculated from both even and odd 
eigenfunctions of the electrostatic potential (j) r with even integer r obtained from Mathieu equation 
(JT3D • The lines colored black, red, yellow, green, blue, and cyan represent r =0, 2, 4, 6, 8, and 10, 
respectively. For odd eigenmodes, r = 2,4, 6, 8, 10. 

these continuity requirements. For odd eigenfunctions <ft r we require r ^ 0. For each value 
of r, the parallel wavenumber K z = k z vo/u pe is chosen to maximize the eigenmode growth 
rate for that r. (See Fig. E]^b,c) below). 

Figure [2] shows Ef 9 and Ef 9 at z = and t = for e = 0.2. For E™ 9 and Ef 9 obtained 
from even <fr r , r=0, 2, 4, 6, 8, and 10; and from odd <p r , r = 2,4,6,8, 10. When z = and 
t = 0, Ey 19 oc —d(p r /dy and -Ef 9 oc r . We scale E e z 19 and i?^ 9 so that their maximum 
amplitudes are normalized to unity based on the assumption that all of the electric-field 
perturbations with different wavelength initially have approximately the same amplitude. 
For even <p r , we can see that E^ 9 vanishes at 9 = 0. By contrast, for even <p r , Ef 9 vanishes 
at 9 = for small values of r = 0, 2, 4 (black, red, and yellow lines, respectively) while E e z 19 
peaks for the larger values of r=8 and 10 (green and cyan lines, respectively). Comparing 
electric fields from odd <p r , E^ 9 is similar to E e z %9 from even <p r , but the peak value at 9 = 
of Ey 19 (r=10) is much lower. E e z %9 obtained from odd r vanishes at 9 = for all r. These 
properties can produce specific features in the 2D electric fields that distinguish the electric 
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FIG. 3. Theoretical growth rate in r-K z space in panel (a) and frequency in panel (d) for case 
e = 0; growth rate in panel (b) and frequency in panel (e) for even <p r for e = 0.2 while growth rate 
in panel (c) and frequency in (f) are for odd (j) r for e = 0.2. respectively. The peak of the growth 
rate in (a) is around (r, K z ) ~ (0, 1), the peak of the growth rate in (b) is around (r, K z ) ~ (10, 0.8) 
and the peak of the growth rate in (c) is around (r,K z ) ~ (2, 1.2). The symbol "+" in panels (e) 
indicate the position of the peak growth rate for e = 0.2 in r-K z space. Panels (g), (h) and (i) 
show ID growth rates vs. K z for cases e = and e = 0.2. The arrows in (g), (h) and (i) indicate 
the directions in which r increases. 

fields in the sheared case from those produced by the instability driven by uniform electron 
drift. For e = 0, the initial electric-field eigenmodes at z = are proportional to either 
cos(r#) or sin(r#) with even integer r. 

We obtained the linear frequency and growth rate of the unstable eigenmodes for both 
even and odd <p r for each pair (r, K z ) for both e = and e = 0.2. We show the results in 
Fig. [31 Panels (a), (b) and (c) show the growth rate (imaginary part of the frequencies) 
in r-K z space, and panels (d), (e) and (f) show the corresponding real wave frequencies. 
Panels (g), (h) and (i) are ID growth rates vs. K z for some specific r values for cases e = 
and e = 0.2, respectively. The arrows in (g), (h) and (i) indicate the direction of increasing 
r. 
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The unstable eigenmodes for the uniform beam in panel (a) (e = 0) are concentrated in 
two distinct areas, representing two different instabilities, one strong and the other weak. 
The stronger instability is the parallel mode at K y ~ 0, K z ~ 1 with peak growth rate 
~ 0.05cu pe , comparable to the cold plasma limit of the growth rate for the fast-growing 
mode of the Buneman instability ^ ma x ~ ^^ pe (i^s-) 1 / 3 ; where tu pe is the electron plasma 
frequency. The growth rate peaks around K z ~ 1, the same as in the uniform cold plasma 
beam limit where v = k z /u pe , and the frequency is oj ~ 0.03o;p e [panel (d)], which is close 
to the cold-plasma limit u r ~ i ^r(^f~) 1/ ' 3 - The peak growth rate of the weaker lower- hybrid 
instability is located near r = 14 and K z = 0.03, with a growth rate of ~ 0.01. Since r = 14 
corresponds to K y = k y vo/co pe ~ 1, we have k z j \Jk1 + ky ~ k z /k y ~ 0.03 ~ a/ m e /mi, which 
corresponds to the ratio of the parallel to the perpendicular component of the wave vector 
of the lower hybrid instability. The corresponding frequency is ~ 0.01w pe ~ u pi [panel (d)], 
also around the cold plasma limit for the lower hybrid instability ui h = u pi /(l + u^/ill) 1 ^ 2 , 
where u P i is the ion plasma frequency and Q e is the electron cyclotron frequency (effectively 
infinite in our model). The growth rate of the Buneman instability decreases steeply with 
K z when it passes the peak, but in panels (g), (h) and (i), we see that the ID growth rate 
curves for r = 6, 8, 10, 12 exhibit a plateau at very small K z < 0.1 in all three cases; these 
plateaus are caused by the presence of the lower hybrid instability. 

Lower hybrid waves are oblique and can interact with both ions and electrons. For the 
ordering f2j <C ulh <C fi e , an approximation in which the electrons are treated as magnetized 
and the ions as unmagnetized is often justified. Under this approximation electrons resonate 
with waves through their motion parallel to the magnetic field, satisfying ou/k\\ ~ Vd- At the 
same time, ions resonate with waves through their motion both parallel and perpendicular 
to the magnetic field satisfying ou/k ~ Vi. The waves resonating with drifting electrons (in 
the ion rest frame) and ions can lead to the growth of both Buneman and lower hybrid 
instabilities, with the former dominated by ion motion parallel to B and the latter by 
ion motion perpendicular to B. Oblique lower hybrid waves therefore transfer the parallel 
momentum of electrons predominantly to the perpendicular momentum of ions. Thus the 
lower hybrid instability can be thought as the oblique limit of the Buneman instability. In 
our model, electrons are assumed to be infinitely magnetized while ions are assumed to be 
unmagnetized and cold, which is the extreme limit of this model. The transfer of momentum 
does not increase the temperature of ions in our approximation. 
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When comparing Fig. |3]^b,e) and (c,f) for the sheared beam s = 0.2, to Fig. |3{a,d) for 
uniform beam, it is obvious that the velocity shear has changed the location of the most 
unstable eigenmodes for even potential eigenf unctions 4> r [panel (b)] to (r,K z ) m (10,0.8), 
with a growth rate close to 0.05u pe . A number of eigenmodes with r < 14 have a comparable 
growth rate. The corresponding frequency shown in panel (e) is ~ 0.04u pe (with the fastest- 
growing mode indicated by the white cross), which is slightly higher than without shear. The 
lower hybrid instability, however, is not affected by the velocity shear, which might be due 
to our assumption of cold ions. The unstable eigenmodes for odd potential eigenfunctions 
4> r in panel (c) with r < 14 have nearly equal growth-rate maxima, with that of the r = 2 
mode being slightly higher than for other unstable modes, but lower than the growth rate 
of the fastest-growing (r = 10) eigenmode for even <p r . Thus we can conclude that the most 
unstable mode is still associated with the Buneman instability. The velocity shear enhances 
the coupling between oblique waves and electrons. The parallel phase speeds for both cases 
are around 0.03-Uo, matching the cold plasma limit. 

III. RESULTS OF VLASOV SIMULATIONS AND COMPARISONS WITH 
THEORY 

We have carried out two 2D Vlasov simulations with two different velocity shears: e = 
and e = 0.2. The simulation box size is L y x L z = 256A e x 800A e . The velocity range is 
[—15, 15]f te . The grids in phase space are (L z , v ez , v iz , L y , v ey , v iy ) = (512, 128, 128, 64, 1, 64). 
The time step is about 0.626cj~ e 1 and The total simulation time is u pe t ~ 1004. We initialize 
the simulations with Maxwellian distribution functions for both electrons and ions, and work 
in the initial ion rest frame. The electrons are given a mean drift vq = 5vt e with a weak 
cosine velocity shear (for e ^ 0) as shown in Eq. (JTJ) . The ion-electron mass ratio is 1836. We 
assume that the magnetic field along z is infinite for electrons, so the electron distribution 
function is restricted to a three-dimensional (z-y-v z ) phase space. The ions are assumed 
to be unmagnetized. The electron-ion temperature ratio is 5 so that the ions are relatively 
cold. 

In §11 Bl we have shown that both Buneman and lower hybrid instabilities can co-exist 
at the same time, but the lower hybrid instability is much weaker. Therefore in the Vlasov 
simulations we do not expect to observe the lower hybrid instability during the linear growth 

11 




FIG. 4. Panels (a) and (b) are, respectively, power spectra of E y and E z in logarithmic scale in 
2D k space at u pe t = 552 for the e = Vlasov simulation; Panels (c) and (d) are the same as (a) 
and (b) for e = 0.2. [Note: some labels are cut off in the embedded PDF version of this figure.] 

phase. In the case of £ — 0.2, the fastest-growing unstable Buneman mode is the r = 10 
potential eigenmode for even <p r . This mode is dominant in the simulation while other 
unstable eigenmodes for r = 0,2,4,6,8 for even <p r compete with the unstable eigenmodes 
for odd (f) r . To verify this we show the power spectra of the electric fields |-E y (k)| 2 and 
|i£ z (k)| 2 in logarithmic scale at u pe t = 552 in both simulations in Fig. HI 

In the case of e = [Fig. 0] (a, b)], the dominant fc-space mode is parallel to z (i.e., 
k y w 0) with K z ~ 1.1. A Fourier transform from the time domain over an interval containing 
oj pe t = 552 reveals the frequency to be about 0.03co> pe . Both the wavevector and frequency are 
consistent with the theoretical predictions for the Buneman instability with uniform electron 
velocity drift. The 2D power spectra [^(k)) 2 and |.E z (k)| 2 for e = 0.2 are shown in panels 
(c) and (d). We see that the strongest unstable modes for E z are now near K z = 0.9, which 
is consistent with the fastest-growing unstable r = 10 eigenmode for even <ft r at K z ~ 0.8 
in Fig. |3]d, and span a large range in K y . The Fourier transform from the time domain 
shows the frequency increases slightly to 0.04o;p e . In this case r has a more indirect relation 
to K y . Here we will use the mapping of eigenmodes from r to K y discussed in Appendix 
IA1 For the globally fastest-growing eigenmode for even eigenfunctions <p r at r = 10 and 
K z = 0.8, the spectrum of E z (K y , K z ) peaks near K y = 0.2. The neighboring E z even 
eigenmodes at r =4, 6, 8, and 12, which have almost as high a growth rate, exhibit multiple 
maxima in their K y spectra, with peaks both at low K y (< 0.2) and high K y (> 0.5) This 
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FIG. 5. The growth in time of the spatially averaged electric field (|-E Z |) in Vlasov simulations for 
e =0 (black solid curve) and e = 0.2 (dashed red curve). 

behavior can be compared with that of corresponding synthetic K z -K y spectra in Figj9] and 
[10] of Appendix A. It is therefore reasonable to conjecture that the strongest modes of E z in 
Fig. H]are produced by the higher-r eigenmodes for even <fr r (e.g., r = 8-12) for which E^ 9 
peaks near the center in Fig. Ek, although the width of the E z spectrum in K z from Fig. HJi 
is slightly narrower than would be expected based on the locus of growth-rate maxima (as 
a function of r) in Fig. [3t>. The modes in Fig. [4] with K z > 1.1, by contrast, come from the 
low-r eigenmodes for even <ft r (e.g., r = 0-4), for which E^ tg peaks near the edge in Fig. [2^. 
The growth rate of unstable eigenmodes for odd <p r compete with the lower r eigenmodes 
for even r and K z ~ 0.8-1.4. These modes spread the spectrum of E z . A similar argument 
can be used to explain the power spectrum of E y . The E z and E y spectra for e = 0.2 peak 
at two distinct ranges of K z , which come, respectively, from the two regions where dv^/dy 
vanishes: the center and edge of the box. At these locations Vb is respectively 20% higher 
and 20% lower than for e = 0. Therefore, the dominant K z in each region should be ~20% 
lower and ~20% higher than for e = if we apply the cold plasma relation k z = u r /vb- 

The time evolution of the spatially averaged electric fields in both simulations shown 
in Fig. [5] can be divided into two stages: linear and nonlinear. In the nonlinear stage the 
electric field reaches its peak and saturates via electron trapping. In the case of e = 0.2, 
electric fields become nonlinear slightly later than when e = due to the small difference 
between the growth rates of the Buneman instability in these two cases. An effective growth 
rate 7 can be determined from the simulation fields during the linear phase through the 
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FIG. 6. Panels (a) and (b) are the Vlasov-simulation electric fields E z and E y for e = at 
u pe t = 679. Panels (c) and (d) are the same as (a) and (b) for e = 0.2 and oj pe t = 753. The two 
times correspond to approximately the same late nonlinear stage. 

relation 7At = ln(E/E ). For E ~ 1 and E ~ 10 -13 , u pe At =600 and 700 corresponding 
to 7 ~ 0.05 and 0.043 for e =0 and 0.2, respectively. These values are consistent with the 
results from linear kinetic theory shown in Fig. El 

In Fig. [6l we show the electric fields E z and E y from Vlasov simulations in z-y space. 
Panels (a) and (b) are for e = at u pe t = 679 and panels (c) and (d) are for e = 0.2 at 
u) pe t = 753, when the simulations for both cases are in a similar (late) nonlinear stage. These 
figures show localized and intense structures, indicative of electron trapping. 

A comparison at a slightly earlier time of structures in z-v z phase space is shown in 
Fig. [71 Panel (a) is the electron distribution function y = 128A e for e = at u pe t = 627, 
(b) is at y — 128A e and (c) is at the edge of y for e = 0.2 at u pe t = 691. Again, because 
the instability with e = grows faster than with e = 0.2, the evolution of electric fields and 
electron holes shown in Fig. [7] are both at approximately the same (early) nonlinear stage. 
That the width of electron holes in (c) is smaller than the width of electron holes in (b) is 
consistent with the wavelength of the Buneman instability in sheared beam. It's interesting 
that the electron holes at the edge seems less regular than the electron holes at the center 
due to the fact that the electric field E z is weaker than the electric field E z at the center of 
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FIG. 7. Panels (a), (b) and (c) are the Vlasov-simulation electron distribution functions in z-v z 
phase space, (a) is at y = 128A e for e = at uj pe t =627; (b) is at the center and y = 128A e and 
(c) is at the edge of the y for e = 0.2 at time ui pe t =691, respectively. The electron holes indicate 
that both simulations are at approximately the same early nonlinear stage. 

To further demonstrate the eigenmode structure of the Buneman instability when e = 0.2, 
we compare 2D electric fields E y , and E z based on the Mathieu-equation analysis to the 
electric fields in the linear stage of the corresponding Vlasov simulation at u pe t = 552. To do 
this, we construct a wavepacket of eigenfunctions with similar growth rates as a superposition 
of eigenfunctions with different weight and phase for r = 0, 2, 4, 6, 8, 10, and 12 for even 
potential eigenfunctions r and r = 2,4, 6, 8, 10, 12 for odd potential eigenfunctions cj) T . From 
the theoretical growth rate, we see that the features of the electric fields come mainly from 
the eigenmodes at r = 10 for even potential eigenfunction r and r = 2 for odd potential 
eigenfunction <f) r . The r = 10 even eigenmode is dominated by small K y (long- wavelength) 
behavior near the center (y = 128A e ) while the r = 2 odd eigenmode is dominated by large K y 
(short-wavelength) behavior near the edges. We draw the 2D E z structure synthesized from 
an eigenmode superposition in Fig. 0(a) and the corresponding E y in (b). Not surprisingly, 
we see that E z peaks at the center (y/X e ~ 128) where the drift is a maximum, while E y 
vanishes there. Compared to the electric fields at u pe t = 552 in the linear stage of the 
Vlasov simulation shown in Fig. [S] (c, d), the superposition of theoretical electric fields 
reproduces the main features of the simulation electric fields, but not necessarily all of the 
subtle details, which are controlled by the exact amplitudes and phases of the different 
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FIG. 8. Panels (a) and (b) show E z and E y , respectively, in z—y space as synthesized from a 
superposition of theoretically determined eigenmodes for e = 0.2; Panels (c) and (d) show E z and 
E y , respectively, from the e = 0.2 Vlasov simulation at u pe t = 552, which is during the linear stage. 



contributing eigenmodes. We show the power spectra of the theoretical electric fields in [TT] 
in Appendix. These features of the electric fields persist into the nonlinear stage, even after 
electron holes form (e.g., at u pe t ~ 700, when the simulation has evolved into the nonlinear 
stage). Even though the electric fields shown in Fig. |6]^c, d) become localized, we still see 
the y-dependence discussed above. On the other hand, E z in panel (a) and E y in panel (b) 
do not show any preferred value of y, as expected for a uniform (i.e., unsheared) electron 
beam. 



IV. CONCLUSION 

In this paper we explored the impact of velocity shear on the Buneman and lower hy- 
brid instabilities. We have studied the unstable modes for electron drifts modulated by a 
weak cosine velocity shear. In this case the dispersion relation can be approximated by the 
well-known Mathieu equation, where we used a kappa-function for the initial distribution 
of infinitely magnetized electrons (with cold unmagnetized ions), and we use a Taylor ex- 
pansion for the velocity shear. These approximations have been validated by comparing the 
analytical results with the results of Vlasov simulations. 

The effect of velocity shear is expected to strengthen the interactions between oblique 
waves and particles because shear causes the fastest-growing Buneman mode to change from 
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a parallel plane wave to an eigenmode containing oblique as well as parallel Fourier modes. 
The shear also distributes the unstable Buneman eigenmodes over a wider range of eigenval- 
ues, with nearly the same slightly lower growth rate than in the shear-free case. We conclude 
that the growth rate, wavelength, and orientation of the fastest growing Buneman mode are 
controlled by the presence (and amplitude) of shear in the electron velocity drift. Therefore, 
velocity shear can lead to different growth rates and wavelengths of unstable modes in the 
direction of the velocity gradient. The Vlasov simulation show that the electron holes form 
at both the center and the edge of the simulation box. At the center, electron trapping 
exhibits a longer wavelength consistent with the fact that the beam is at its fastest; at the 
edge, electron trapping exhibits a shorter wavelength, consistent with the correspondingly 
slower beam. 

The velocity shear does not affect the growth rate of the weak lower hybrid instability, 
in which the interactions between lower hybrid waves and ions can transfer electron parallel 
momentum to the perpendicular motion of the ions. This momentum transfer can occur 
without necessarily affecting the temperature of the ions. 

The advantage of our method is that it provides eigenfunctions for every unstable mode so 
that we can investigate both the spectra and 2D spatial structure of the electric fields. These 
results can help us to interpret the physical content of simulations, experiments and satellite 
data. Our method has assumed that the shear is weak. However, stronger velocity shear may 
reveal new physical regimes, including shear-driven as well as shear-modified instabilities. 
To address this issue, we are undertaking a study of how far the Mathieu-equation analysis 
can be extended. 
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Appendix A: Mapping Mathieu equation eigenfunctions into Fourier space 

We build a map between r and k y by projecting the complex eigenmodes of E y (y, k z ) and 
E z (y, k z ) at t=0 into the z-y plane, 
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FIG. 9. The correspondence between K y and even integer r. Panel (a) shows E y and panel (b) 
shows E z for e = 0.2. The white stars mark the linear correspondence K y = m/L y for e = 0. 



Ef 9 = -Re(ik z (cos(k z z) + i sm(k z z))(Re((f) r ) + ilm((j> r ))) 

= k z cos(k z z)Im(<j) r ) + sxn(k z z)Re(<p r ) , (Al) 

7T d 

E v 9 = --^Re(ik z (cos(k z z) + ism(k z z)) — (Re((f) r ) + ilm((fi r ))) 

Ly do 

= ^-{K cos(k z z)Im(d(j) r /d6) + sm{k z z)Re{d<t) r /d6)). (A2) 

Ly 

We then perform a 2D Fourier transformation of E^ 9 and E^ 9 in (K y , K z ) = vo/u pe (k y , k z ) 
space. The results of Fourier transformation should be the same for even and odd eigen- 
functions. In this appendix, we use even eigenf unctions <p r to analyze the relation between 
r and K y . The results show that for e = 0.2 the corresponding power spectrum is no longer 
concentrated on one specific K y but is instead spread over a range of K y , as shown in Fig. |9j 
We take the K y corresponding to the maximum power intensity. On the other hand the 
correspondence between r and K y is different for E y and E z . Fig [TOl shows this map between 
K y and r. Black circle is for E z . The E z eigenmodes r = 8, 10, 12 have K z ~ 0.8 while 
r = 10 corresponds to K y ~ 0.2, r = 6, 8, 12 correspond to K y ~ 0.6, 0, 0.5 respectively, 
r = eigenmode corresponds to K y ~ 0.1, K z ~ 1.1 and others in between K z ~ 0.8 — 1.1 
and K y ~ 0.1 — 0.6. The E y eigenmodes for r = 0,2,4,6,8,10 have a relatively linear 
correlation with K y , but for r = 12, K y decreases from 0.9 to 0.5. 

To provide a fuller picture of the spectra from theory, we show in Fig. [11] the power 
spectra of the superposition of theoretical eigenmode electric fields E y and E z previously 
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FIG. 11. The power spectra for theoretical electric fields E y and E z shown in Fig. [8] (a, b). Panel 
(a) is for E y and panel (b) is for E z . 

shown in Fig. [HJ We reiterate that the 2D theoretical E y and E z only reproduce the main 
features of the electric fields and thus the their power spectra shown in Fig. [TT] can only 
approximate the realistic power spectra from the simulation as shown in Fig. HI However, 
Figs. H] and [11] do both show the same two distinct ranges of K z . Again, the differences are 
controlled by the exact amplitudes and phases of the contributing eigenmodes, which we 
can not determine theoretically. 
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